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ABSTRACT 


The need of telecommunications continue to increase and new methods of utilizing 
more of the potential bandwidth through optical fiber become more popular. So the 
demands placed on light wave systems have dramatically increased, spurring 
development of light wave. components. So the demand of modeling of optoelectronics 
devices has increased the importance of numerical methods in sharp manner. Beam 
Propagation Method is an very efficient tool to analyze the propagation of wave through 

p, 

the complex structure. 

In this Thesis Finite Difference Beam Propagation Method is used to analysis the slab 
wave guide. Slowly varying envelope Approximation is used which help the method to 
choose the relatively bigger step size, lesser memory space, lesser computational time as 
compared to other numerical technique. Transparent Boundary Condition is also used 
because of limited computational window. Analysis of perturbed Slab wave guide 
structure is done. Effect of cladding removal and effect of different Refractive Index at 
the substrate in slab wave guide structure is also analyzed. 
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Chapter 1 

INTRODUCTION 

1.1 Introduction 

The use and requirement of telecommunication continues to increase and new 
methods, utilizing most of the potential bandwidth through optical fiber is becoming 
more popular. The demands placed on light wave systems have dramatically 
increased, spurring development of light wave components such as fibers, lasers, 
detectors, modulators, switches, and wavelength multiplexing/demultiplexing devices. 

The analysis of theses devices is useful for the development of optical 
communication. Nowadays many types of analysis of such components exist: Exact 
Analytical, Numerical and Semi-Analytical methods. Exact analytical methods are 
heavily restricted to only very simple slab devices, limiting their use in the analysis of 
complex 3D structures. Consequently, in the quest for analysis, more attention has 
been focused on Numerical and approximate Semi-Analytical methods. 

1.2 Beam Propagation Method (BPM) 

Among the different numerical techniques developed for modeling guided-wave 
propagation in integrated optics and fiber optics, the most popular is the beam 
propagation method (BPM)[1][2]. Unlike the finite-difference time-domain (FDTD) 
technique and finite element Method (FEM), which simulates the device's 
performance through a rigorous numerical solution of Maxwell's equations, the BPM 
takes advantage of some simplifications common to many problems. 
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The simplifications offered by the BPM reduce the computational complexity of most 
problems, speeding up the computational process for many problems, and in some 
cases making otherwise impossible computations realistically solvable. By 
eliminating the second derivative term in z, the problem of solving Helmholtz 
equation is reduced from a second-order boundary value problem requiring iteration 
or eigenvalue analysis to a first-order initial-value problem that can be solved by 
simple "integration" along the propagation direction z. 

Numerical methods such as FTDT or FEM sometimes require large amounts of 
memory and computational resources. Consequently, in certain situations, it is 
extremely difficult to apply numerical methods, especially where the structure is large 
and complex (such as the array waveguide structures). In these cases, a Semi- 
Analytical method [3], such as the Effective Index (El) method [4], Spectral Index 
(SI) method [5], or Free Space Radiation Mode (FSRM) method [6] may be chosen 
for analysis, depending on the geometry. Each method uses approximations to reduce 
the complexity of the device, which can result in a fast application suitable for design 
purposes, with the additional advantages of much lower requirements in CPU time 
and resources. 

In this thesis research numerical method. Finite Difference Beam Propagation 
(FDBPM), is used for the analysis and design of evanescent wave based guided wave 
optical devices for the measurement of refractive index of the liquids. Propagation of 
monochromatic wave in slab waveguide structure is studied. Here 2D Waveguide 
structure is used for the analysis of the effect of refractive index. The Structure of 
waveguide is as shown in figure 1.1. 



, C 
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Window 


Figure 1.1 Structure of Slab waveguide 

In this Thesis research the output power is measured with change in refractive index 
of the Surrounding medium as well as change in different amount of the cut in the 
cladding. Effect of concentration of sucrose solution on the output power is also 
analyzed. Effect of perturbation and strips on core is also studied. 
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1.3 Thesis Organization 

Chapter 2 describes the theory of BPM and FDBPM. In this chapter mathematical 
model of finite difference method, Transparent Boundary Condition will be discussed 
and applied into numerical method. 

Chapter 3 deals with the accuracy and validation of simulation program. 

Chapter 4 discusses the effect of different percentage of cross sectional cut in the 
slab waveguide structure. In this chapter we also see the effect of change of refractive 
index at the surrounding medium of slab is also discussed. 

Chapter 5 concludes the thesis with the major results and gives suggestions for 


further work in this area 



Chapter 2 

FINITE DIFFERENCE BEAM PROPAGATION METHOD 
THEROY AND ITS IMPLIMENTATION 
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In this chapter, the basic concepts of Finite Difference Beam Propagation 
Method (FD-BPM) [7] is discussed and its algoritlim is derived. This derivation is 
based on the Maxwell’ Wave Equations [8][9], vector wave equation, and its the 
Scalar form is derived. The finite difference method is used with Crank-Nicholson 
(CN) propagation scheme [10], known for its accuracy and stability, to solve the 
scalar form of wave equation. Transparent Boundary condition [11] is applied to 
reduce the effect of the simulation space limits on the solution, such that it can be 
used for a real time application. Methods of implementation into a computer program 
are explored with the aim of making a set of reliable and efficient solvers that are 
suitable for the simulation of large components. 

2.1 Finite Difference Beani Propagation Method Principles. 

I 

xt 


Surround Medium 


Cladding 

Z 

Core 

Cladding 


Surround Medium 


Figure 2.1 Slab Waveguide structure used for thesis 
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The Finite Difference - Beam Propagation Method (FD-BPM) [7] is a popular 
numerical algorithms used to determine light propagation in the physical structures 
like slab waveguide (Figure 2.1.)[12][13]. It is a used to solve the MaxwelFs Wave 
Equations and it can be applied to complex waveguide structures. A good commercial 
example of BPM software package is BeamProp by RSoft Design Group, Inc. [14]. 
The FD-BPM samples both the structure and field in analysis. The structure was 
uniformly sampled in steps of Az in the direction of propagation of the wave and the 
field was also uniformly sampled in steps of Ax in the transverse direction of 
propagation of wave. Figure 2.2 demonstrate the FD-BPM working in the 2D case. 
From a known set of ‘present’ points, forming a line in 2D implementations as shown 
in figure 2.2, the field at the ‘next’ set of points is calculated. Knowledge of the 
structure at all points is required. This process is repeated for each new set of known 
field points, and so the simulation ‘propagates’ in the z-direction. 


Present ! 

Sample points i Sample Points 

X ! X 


X I X 


n, 




Transverse 
Direction (x) 



X I X 
X i X 


Propagation 
Direction (z) 


Figure 2.2 The BPM algorithm calculates the ‘next’ field from the ‘present’ field. 
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2.2 Derivation of the FD-BPM Algorithm 
2.2.1 Maxwell’s Equations 

The following Maxwell’s wave equation are solved by using Finite Difference — 


Beam Propagation Method (FD-BPM) [7]. 

In a source free region, and assuming a 

periodic time variation, they are: 


Ax£= 

(2.1) 

Ax H=jfj.coE 

(2.2) 

A»sE = 0 

(2.3) 

=0 

(2.4) 


where e=^ereo and = and the vector quantity E (V/m) is the electric field 

vector and H (A/m) is the magnetic field vector and a periodic time variation e"” has 
been assumed. The quantities s and p define the electromagnetic properties of the 
medium. £o = 8.854 x lO"'^ Fm'^ is the dielectric constant in vacuum and 
//o = 4 ;r xlO'^ Hm'^ is the magnetic permeability in a vacuum. Sr and pr are there 
relative permittivity and permeability of the material. In this work pr has unity value 
since only non-magnetic materials are considered, er had the different value 
depending on the type of different surrounding material used in the analysis. Using Sr 
and pr the refractive index n, of the medium is defined as n = -^jLirSr . 

2.2.2 The Wave Equation 

The total electromagnetic field that is supported by a waveguide can be expressed in 
terms of only the electric or magnetic field components to produce wave equations. 
Throughout this section our calculation considers only the electric field, rather than 



the magnetic field. This was done by taking the curl of equation (2.1) and substituting 
in (2.2), i.e. 




(2.5) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


For TE Mode f Ex = Ez = Hy =0) 
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Substituting value of Ex, Ez, Hy in equation(2.5), equation(2.6), equation(2.7), 
equation(2.8) and equation(2.9) 


dEy . „ 

- -jcoixoHx 

dz 

(2.11) 

dEv 

— = -jcopioHz 
dx ^ 

(2.12) 

dHx dHz 

“ — jCOSoSrEy 

dz dx 

(2.13) 

Putting the value of from equation(2.1 1) and the value of Hz from equation(2.12) 

in equation(2.13) 


a'£v d^Ey , ^ „ 

dz^ dx^ ^ ■ 

(2.14) 

TM Mode (Hx= Hz = Ey = 0) 


Substituting value of Hx, Hz. Ey in equation(2.5). 

equation(2.6), equation(2.7). 

equation(2.8) and equation(2.9) 


dEx dEv , 

-r ^ = -JCOfloHy 

dz ox 

(2.15) 

dHy . „ 

— J (DSoSrEx 

dz 

(2.16) 

dHy . „ 

= JCOSoSrEz 

dx 

(2.17) 


Putting the value of Ex from equation (2.16) and the value of Ez from equation (2.17) 
in equation (2.15) we gets 

^_|l + A(J_A(^,Ex))+ko^erEx=0 (2.18) 

az' dx Sr dx ^ ’ 

Where ko is wave number in free space, X.o is wavelength in free space. 


(2.19) 
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2.2.3 FD-BPM Formulation 

The uniform sampling discussed previously is used to solve 
second order differential equation (equation 2. 14). 

d^Ey d^Ey , ^ ^ 

7 -H Z- + CO,r£oU„Ey =0 

dz^ dx^ 

Assuming the solution of equation (2.14) is as follows 

Ey(x,y,z)=O(x,y,z)exp(-y;0z) (2.20) 

Slowly envelope approximation is used to separate Ey into slowly varying envelope 
function <E)(x, 7 ,z)and fast oscillating phase term ex\>{-jPz). /3 is propagation 
constant in z direction. [3 and koare related as follows. 

fi = n<!ti"x.ko ( 2 . 21 ) 

where na/r is refractive index. 

Putting the value of Ey in equation (2.14) we get 


5^0 


- 2yy5z— -/?^(D+-^+ko^frO=0 


ao 

dz 


dz^ dz ^ ^ ■ obc' 

As 0(jc,>^,z) is a slowly varying envelope function in z direction 


ao a'o 

— » — r 
dz dz' 


( 2 . 22 ) 


So 


dz^ 


= 0 


(2.23) 

(2.24) 


Using this equation (2.22) becomes 


a^o ao -) ^ 

£21 -2jpz- /?^ 0+ko^£rfl5=0 

dx ' az 


(2.25) 



II 


2.3 Finite Difference Method:- 

Finite difference method is an idea of approximating the derivatives of a continuous 
function by dividing up the function into segments and looking at how the slope (first 
derivative) varies between segments. Consider a function u{x), divided along the x- 
axis into segments of width Ax. Consider one grid point (segment boundary) near the 
center at grid point at position x given by xj. The point to the left at x-Ax is grid point 
Xj.\ and the point to the right is at x+Ax is at grid point xj+\. We can extract 
information about the first and second derivative of the function u{x) from the 
information and grid points by using the following technique. 

2.3.1 Taylor Series Expansion 

The Taylor series expansion relates the value of the function u{x) to values at adjacent 
points (x+Ax or x-Ax) and the derivatives of m ( x ). For example; 


u{x + Ax) = w(x) + Ax — + 


du Ax^ d^u Ax^ d^u 


and 


dx 2 dx^ 6 dx^ 


. . . , ^ ^ du Ax^ d^u Ax^ d^u 

u{x-Ax) = w(x)- Ax H r r + 


(2.26) 


(2.27) 


5x 2 dx^ 6 9x’ 

Thses equations are used to determine the first order and second order derivative in 
the following section. 


2.3.2 First Derivatives ( — ) 

dx 

Forward Difference 

If we keep only the first two terms of the expansion for m(x+Ax), then 

du _ m(x + Ax) - m ( x ) 
dx Ax 


(2.28) 
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This is a first order accurate solution {0 (Ax)). It is called the forward difference since 
it finds the derivative in the positive (forward) direction of the profile from the point 
atx. 

Backward Difference 

Considering only the first two terms of the expansion for u(x-Ax), in equation (2.27) 

du _ u(x)-u(x + Ax) 
dx Ax 

This is a first order accurate solution (0(Ax)). It is called the backward difference 
since it finds the derivative in the negative (backward) direction of the profile from 
the point at x. 

Centered Difference 

To get a more accurate representation of the first derivative, find the difference 
between w(x+Ax) and u{x-Ax) is combined using equation (2.26) and equation (2.27) 
we get 


du Ax^ d'u du Ax^ d'u 

u{x + Ax) - u{x - Ax) = m(x) + Ax — H - m(x) + Ax- 


dx 2 dx' 
On simplification we get the following equation. 

du w(x + Ax)- u(x- Ax) 


dx 2 dx^ 


(2.30) 


(2.31) 


dx 2 Ax 

This is a second order accurate solution (here the Ax^ terms is included, they just 
happen to cancel out - 6)(Ax^)). It is called a centered difference, because it is centered 
at the position x, but uses information from both adjacent points. 


2.3.3 Second Derivative ( 


d-u 

dx^ 


) 


For the differential equation, we actually need the second derivative. Using the Taylor 
series expansion, the sum of h(x+Ax), form equation (2.26) and, i;(x-Ax) form 


equation (2.27) we get 
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dx 2 dx^ 


dx 2 dx~ 


On simplification we get 


u{x + Ax) + u{x - Ax) = 2u{x) + Ax' 


d^u 

ax' 


Solving for 


d-u 

dx' 


(2.33) 


d'u _ u(x + Ax) + u{x-Ax)-2u{x) 

^2 A '> (2.34) 

OX Ax' 

This is again a centered difference, but it also uses information about u(x) at position 
X. It is also second order accurate( 0(Ax')). 

The accuracy can be improved by decreasing Ax. This is the easiest and fastest way to 
improve accuracy, but it comes with a computational cost and taking more time and 
memory. 


2.3.4 Numerical application of finite difference method 

Again using equation (2.25) 

5<D , , 

2/7?z — = -^-/0'cI)+koVrCl> (2.35) 

dz dx' 

Using discretization in x and z direction we get 

x = pAx (2.36) 

z=lAz (2.37) 

p, 1 are integers. Az and Ax are uniform small steps in the direction of propagation and 
transverse to the direction of propagation respectively. Figure 2.4 shows the 

discretization scheme on the computational window 
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Z, 

Z2 
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Z3 
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Figure 2.4 Discretization in the computational window 


following notations are used for wave function 0(jc,z) and relative 
permittivity 5r(x,z) . 


0(;c,2) = 0'/. 

(2.38) 

^■/( a :,^) = St' ( p) 

(2.39) 


Using this in equation (2.34) we get the discretization in the z direction as follows 


5 0 1 . 0/7 + 1 0/7 


dx^ 


/Sx 


Ax 


AiX 


-) 


(2.40) 


0/7 + I — 0/7 

Aa 


0/7 0/7 - 

Ax 


has centre difference at p+1 /2 

■ has centre difference at p-1/2 

5^0 . 0/7 + 1 — 20/7 + 0/7- K 

dx^ Aa^ 

Equation (2.41 ) has centre difference at point p 

ko^(er - ne//' )0 = ko^(Sr(p) - neff')^p 


(2.41) 


(2.42) 



Equation (2.42) has centre difference at point p 

Putting the values in equation(2.35) form equation(2.40) and equation(2.42) we get 
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dz 


f I' p I p - I, ,2/ 7/ N 

2jp — = { — ) + ko {sr {p)-nerW P (2.43) 


Ax" 


Now discretizing the equation (2.35) in z direction. 

d<i)n o''"' n - o' « 

2jl3^ = ljl5{ ~ \ 

oz Az 


(2.44) 


Discretizing centre of left hand side of equation (2.44) is at /+1/2 while right hand 
side of equation (2.35) has discretizing centre at /. 

Now using Crank Nicholson (CN) scheme to get the solution of equation (2.35) at the 
same discretizing centre with good stability. If we have field profile at z and at z + Az 
then its discretization center will be at z + Az/2. Crank Nicholson (CN) scheme has 
proposed an difference parameter a . Which defines the a amount of field at z and 
(1- a) amount of field at z + Az are added to calculate the field at z + Az/2. For 
good stability a is proposed 0.5. 

So equation (2.41) will be as follows. 


2jl3{ 


Az 2 


Ax" 


) + fco"(£r'^' (/>)-«, /r)0'"';,] + 


Ax 




2Ax' 


Az 2Zkx‘' 2 


2Ax" 


(t)' - 1[ , ] + <t>' p{ {Sv {p)- ne/f ' )] + <I>' /- + 1 [ ,] (2.45 ) 

2Ax‘ Az Ax‘ 2 2Ax' 


Let OCw = 


Ax" 


ax = - 


Ax" 
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ae 


1 




tsz 


- 

6c jc ■ 


Ao" (f (p) - neff ^ )] - = 


cD'p - i[a'w] + + a'. + A:o' (^/ (p) - )] + + i[a^] 

Ziz 


(2.46) 


2.4 Selection of «e^ 


P 


fif) 

Figure 2.5 co-p curve 

In the Slab Wave structure only those mode will propagate for whom The ffl/p is fix. 
During the analysis of slab waveguide structure the we have used the single mode 
structure. Single mode have no cutoff frequency. While second order mode have 
cutoff frequency as shown in figure 2.5. So we have taken the value of rieff 5 percent 
of difference of core and cladding refractive index value less from the core refractive 



index. 
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2.5 Transparent Boundary Condition 

An infinite wide area have no analysis boundary and there fore no reflections are 
coming from the boundaries. But such an area cannot be assumed in actual design, 
therefore, a limited analysis window is used. In actual structures, radiated waves are 
reflected at boundaries and return to the core area, where they interact with the 
propagating fields and degrades the calculation accuracy. In this section TBC is 
discussed which is developed by Hadley[IO]. It is a way to efficiently suppress the 
reflection at the boundary. 


p=0 

Mm m m 

p=l 

p=2 

p=3 

p=M-l p=M 

p=M+l 

m mm M 

w" ■ 

Xo 

X, 

X2 

X3 

Xm-1 

m m m 

Xm+1 


Figure 2.6 Nodes at p=0 and p=lVl+l are outside the analysis window 

Assuming that the analysis window (figure 2.6) contains nodes from p=l to p=M. To 
analysis the effect of reflection at the computational boundaries hypothetical nodes at 
p=0 and p=M+l are assumed outside the analysis window. 

In the following sections Transparent Boundary Condition (TBC) is applied at the left 
hand boundary (p=l) and at the right hand boundary (p=M) of the analysis window 

2.5.1 Left Hand Boundary: 

Considering the left hand boundary in figure 2.5 we incorporate the influence of the 
hypothetical node at p=0(outsidc the analysis window) into the node at p=l (inside 
the analysis window) The wave function for the left traveling wave with the x- 
directed wave number kx is expressed as 
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<I)(x, z) = A(z) exp(Jk.a:) (2.47) 

Assuming that the notation for x coordinate and the fields of the nodes at p=0, 1, 2 is 
xo, X|, X2and <J>o, respectively it is assumed that 

Oo(x, z) = A(z) exp(JkKXo) (2.48) 

<I>i(x,z) = y4(z)exp(y^a:i) (2.49) 

02(x,z) = A(z)exp(jkxX 2 ) (2.50) 

The field O i and <4) 2 are inside the analysis window while O 0 is the hypothetical field 
whose influence is incorporated into the field inside the analysis window. 

Dividing equation (2.50) by equation(2.49) we get 
02 


Oi 


■ exp(J/crAx) 


(2.51) 


Dividing equation (2.49) by equation (2.48) we get 
Oi 


Oo 


= exp(JkxAx) 


(2.52) 


Where Ax = X2-X1 = Xi-Xo 
O2 


Let 


O: 


= exp(7A;vAx:) = 71 


(2.53) 


Now using equation (2.53) and equation (2.52) we get 



71 

Now the value of x directed wave number is obtained from equation (2.53) 


(2.54) 


L = — ln(7i) (2.55) 

JAx 

Since wave travels leftward the real part of the x-directed wave number b should be 
negative. When it is positive, which implies reflection at the left-hand boundary, the 
sign should be changed from plus to minus. 
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Oo = Oi exp(7LAx) (2.56) 

2.5.2 Right Hand Boundary 

Considering the right hand boundary in figure 2.5 we incorporate the influence of the 
hypothetical node at p=M+l (outside the analysis window) into the node at 
p=M(inside the analysis window) 

The wave function for the left traveling wave with the x-directed wave number kx is 
expressed as 

<I>(x,z) = .4(z)exp(-yA;dc) (2.57) 

Assuming that the notation for x coordinate and the fields of the nodes at p=M-l, M, 
M+1 isxM-i, XM.XM+tand <Dm-i,<1>m.^m+i respectively 
It is assumed that 

<D>w - i(jc, z) = A(z) exp(-JkxXM - 1 ) (2.58) 

z) = A(z) exp(-j/cxXM) ( 2 . 59 ) 

(t>w + i(x, z) = A(z) exp(-jkxXM ♦ i) (2.60) 

The field <I>M-iand<I)M are Inside the analysis window while Om+i is a hypothetical 
field whose influence should be incorporated into the field inside the analysis 
window. 

Dividing equation (2.60) by equation(2.59) we get 

= exp(-JkxAx) (2.61) 

Dividing equation (2.59) by equation(2.58) we get 

— ■ =exp(-JkxAx) (2.62) 

Ow - 1 


Where Ax = Xm-xm-i = Xm+i-Xm 
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Let 


(E>M - I 


e\p{-jLAx) = rjM 


(2.63) 


Now using equation (2.61) and equation (2.63) we get 

Om + 1 = (t>w X T]M (2.64) 

Now the value of x directed wave number is obtained from equation (2.62) 


kx = ■ 


1 


jhx 


ln(7M) 


(2.65) 


Since wave travels rightward the real part of the x-directed wave number kx should be 
negative. When it is positive, which implies reflection at the right-hand boundary, the 
sign should be changed from plus to minus. 

O/w + \ = exp(-_/^:vAjc) (2.66) 


2.6 Numerical Application of TBC 

We know the field O^at z and we want to know field <I> at z-i-Az. Recalling 
Equation (2.46) 

- i-a'*' .] + O'"' p[^ -a'^'x-ko- (e/"' (p) - nar )] - O'"' . i[a'"'c.] = 
iSz 

^'p - l[a'w] + O' p[^ -h a'x + ko^ {Sr' ip) - He/ )] + O';, + l[a'e] 

Az 

where unknown field at z-i- Az are O ^"'p.i,0 ^^'p,0 ^^'p+i and known field at z are 
O ^p.i,0 ^p,0 ^p+i let us assume known coefficients A(p), B(p), C(p), D(p), E(p), 
F(p) to be 


+ 

1 

II 

< 

(2.67) 

B(p) = {'^-^^ a'^'x ko\s/^\p)-nc/)} 

Az 

(2.68) 

C(p) = 

(2.69) 
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D(p) = a\r (2.70) 

E(p) = a'. +^ + k<f {e,'{p) - nc/ ) (2.7 1 ) 

tsz 

F{p) = a'. (2.72) 

Now using these equation, equation (2.45) will be 
0'"'/,-i[A(p)] + cI)'^';,[B(p)] + cD'"'p+i[C(p)] = <D'p-iD(p) + cI)';,E(p) + <I>'p + iF(/7) (2.73) 
Where p represents the lateral position of the node 


2.6.1 Left Hand Boundary (p=l) 

The field <[> i of the node on the left hand boundary 
is influenced by the field Oo of the hypothetical node outside the analysis window. 


Where 


d)o = (Di*YA 



1 


(2.74) 

(2.75) 


Parameter Yi. is determined by the known field at x(i.e. /), Yl is known so 
equation(2.73) can be reduced to 

O'" i[B'( 1 )] + O'" 2[C(1 )] = o' i[E’(l )] + O' 2 [F( 1 )] (2.76) 


where 


B’( 1 ) = - a'"x - U (£r'" (1) - «.//')} 

tsz 


/+! 


C(l) = -a 

E'( 1 )=«'>.■ • yi + a! X + ^ (1) - rui ) 

Az 


(2.77) 

(2.78) 

(2.79) 


F(l) = a'. 


(2.80) 
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2.6.2 Right Hand Boundary (p=M) 

The field d>Mof the node on the left hand boundary 
is influenced by the field Om+i of the hypothetical node outside the analysis window. 


Ow + 1 = Ow * T« 


(2.81) 


Where Yk = t]m = — ■. 

iW / C> A/ - I 


(2.82) 


Parameter Yr is determined by the known field at x(i.e. /), Yr is known so 
equation(2.73) can be reduced to 


- i[A(M)] + 0'^V[B'(M)] = O'm - :[D(1V1)] + a)'M[E'(M)] (2.83) 


Where ^(M) = 


(2.84) 


B’(M) = - a'^'. - (M) - «c-//^)} 

Az 


(2.85) 


D{M) = a'« 


( 2 . 86 ) 


EXM) = Y/ia'c. + a'. + ^ + ko^ (sr'iM) - nar ^ ) 
Az 


(2.87) 


Rewriting equation (2.73) with constants we get the matrix form 


■fi'tl) 

C(l) 

0 

0 

0 




'£’(1) 

£(!) 

0 

0 

0 


A(2) 

5(2) 

C{2) 

0 

0 


(D2'^‘ 


0(2) 

£(2) 

£(2) 

0 

0 


0 

-4(3) 

5(3) 


0 

X 


= 

0 

0(3) 

0(3) 


0 

X 

0 

0 



C(A/-1) 




0 

0 



F[M-\) 


0 

0 

0 

A{M) 

B'{M) 




_ 0 

0 

0 

D(M) 

E\M) 




'b\\) 

C(l) 

0 

0 

0 


-4(2) 

5(2) 

0(2) 

0 

0 

Let BW = 

0 

-4(3) 

5(3) 


0 


0 

0 



C(A/-I) 


0 

0 

0 

A{M) 

B\M) 
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and FW = 


£'(l) 

0 ( 2 ) 

0 

0 

0 


0 ( 1 ) 

0 ( 2 ) 

0 ( 3 ) 

0 

0 


0 

0(2) 

0 ( 3 ) 

0 


0 0 
0 0 
0 

D(M) E\M) 


Then field at the next step z+ Az will be 




'<D/ ‘ 

cDa'^' 

= [BW\\fWx 

CD 2 ' 



0m' 


( 2 . 88 ) 


From the equation (2.88) by the knowledge of matrixes Bff , FW and field O' we 
can calculate the field at the next step. This process can be used for certain distance z. 



Chapter 3 

ACCURACY OF THE SOFTWARE 
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Accuracy is primary and essential condition for any software to be valid for further 
use. To prove the accuracy of this software, one comparison is made with Fresnel 
equation. In this chapter first free space propagation of Gaussian wave is described 
and its Full width at half maxima (FWHM) is measured. As FWHM can completely 
describe the shape of the Gaussian profile so FWHM is used to measure the accuracy 
of the software. Then applying the same excitation as an input to the simulation 
program and FWHM is measured at same distance. At the end error is calculated 
between these two FWHM. It is repeated for different spacing (Az). Among all these 
error, the minimum one is considered for the further calculation of the different 
waveguide structures. 

Initial conditions considered to make comparison is as follows:- 

• Computational window size = 5 1 pm 

• Wavelength = 1 .00 pm 

• No of points in the x-axis = 500. 

• Initial excitation = Gaussian Profile 

• Beam waist of initial Gaussian profile = 3 pm 

• Amplitude of Initial Gaussian Profile = 1 Volt 

• Refractive index of the medium=I .00 

3.1 Laser Beam Propagation 

In general, laser-beam propagation can be approximated by assuming that the laser 
beam has an ideal Gaussian intensity profile, corresponding to the theoretical TEM^^p 
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mode. The TEMp^mode output of a laser is a spherical Gaussian beam. The wavefront 
normals (rays) of a plane wave are parallel to the direction of the wave so that there is 
no angular spread, but the energy extends spatially over the entire space. The 
spherical wave, on the other hand, originates from a single point, but its wavefront 
normals (rays) diverge in all directions. Waves with wavefront normals making small 
angles with the z axis are called paraxial waves. They must satisfy the paraxial 
Helmholtz equation. A solution of this equation is a wave called the Gaussian beam. 
The intensity distribution in any transverse plane is a Gaussian function centered 
about the beam axis. The width of this function is minimum at the beam waist and 
grows gradually in both directions. The wavefronts are approximately planar near the 
beam waist, but they gradually curve and become approximately spherical far from 
the waist. 


The Gaussian Beam 

A paraxial wave is a plane wave exp(-yAz) (with wavenumber k = 2nlX and 
wavelength 2 ) modulated by a complex envelope A(r) that is a slowly varying 
function of position. The complex amplitude is [16] 

U(r) = A(r)exp(-Jkz) (3.1) 

For the complex amplitude WfrJ to satisfy the Helmholtz equation, 
V^U + k^U = 0,the complex envelope AfrJ must satisfy the paraxial Helmholtz 
equation [16] 

AiA-J2/t~ = 0 (3.2) 

dz 


Where A]. = 


dx^ dy2 


is the transverse part of the Laplacian operator. One solution 


to the paraxial Helmholtz equation is paraboloidal wave for which 
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A o 

v4(r) = — exp(-/A— ), where = x" + 3.3) 

z 2z 

The paraboloidal wave is the paraxial approximation of the spherical wave 
f/(r) = {AI r)exp(-jkr) when x and^* are much smaller than z .Another solution of 
the paraxial Helmholtz equation provides the Gaussian beam. It is obtained from the 
paraboloidal wave by use of a simple transformation. Since the complex envelope of 
the paraboloidal wave equation (3.3) is a solution of the paraxial Helmholtz equation 
(3.2), a shifted version of it, with z-^ replacing z where is a constant, and the 
solution is 


A{r) = 


A. 

9(2) 


exp 


-7^ 


2q{z) 


q{z) = z-^ 


(3.4) 


This provides a paraboloidal wave centered about the point z = ^ instead of z = 0. 
When ^ is purely imaginary(^ = -jzo , where zo is real,) equation (3.4) gives rise to 
the complex envelope of the Gaussian beam as follows 


.4(?') = Aexp 

q{z) 


-Jk 


2q{z) 


,q{z) = z+j\ 


(3.5) 


The parameter zo is known as the Raleigh range. 

We can write the complex function l/q(z) = 7/(z + jz) in terms of its real and 
imaginary parts as 

1 \ . X 


J- 


(3.6) 


q{z) R(z) ■' nW^{z) 

W(z) and R(z) are beam width and wave front radius of curvature, respectively. An 
expression for the complex amplitude U(r) of the Gaussian beam is obtained by 
Substituting equation. (3.6) into equation.(3.5) and using equation.(3.1) as 


V{r) = A, 


fV(z) 


exp 


r 1 


W\z)_ 

exp 


Jkz-jk 


2R{z) 


+ y^(^) 


(3.7) 



The beam parameters can be evaluated by following equations 


W{z) = W, 


R{z) = z\ 


^(z) = tan' 


■ 

f \ 

2“ 

1 + 

z 


— 






1/2 


■ 

( r \ 

2 “ 

1 + 

fiL 



1 ^ J 





V^oy 


W. 


f— 1 


1/2 


(3.8) 


(3.9) 


(3.10) 


(3.11) 


Now consider a beam wave Uo{r) propagating in the z direction in free space .At the 
transmitting aperture (z=0), the wave has a Gaussian amplitude distribution with beam 
size Wo and a curved phase front with radius of curvature i?o.Then 


t/o(0,p) = exp 


1 . k 


(3.12) 


In . 


where k = — is wave number, A is wavelength, p is the radial vector from the z axis. 


--^kap^ 


We can also write Equation. (3.12) as [15] 

C/o(0,/7) = exp 

. . T . 1 

where a = a. + /a, = ;-+ / — 

at an arbitrary point (/C>,z) the beam in free space is given by 

1 


(3.13) 


^o(A^) = 


(1 + jaz) 


exp 


jkz- 


ka\ 


(l + yoz) 


(3.14) 


this is valid within a distance z [15] such that 
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which covers almost all the laser propagation distances used in practice. The effect on 
the beam parameters as beam advances is described below. 


Intensity 

The optical intensity I{r) = \U{r)\' is a function of the axial and radial distances z 


and p = {x' + y 



I{p.^) = h 


‘ K 1 

2 

2p^ ■ 

_lT(z)_ 

exp 

fV^(z)_ 


(3.16) 


Where /() = . At each value of z the intensity is a Gaussian function of the radial 

distance p . The Gaussian function has its peak at /? = 0 (on axis) and drops 
monotonically with increasingp. The width JV(z) of the Gaussian distribution 
increases with the axial distance z as illustrated in Fig.3. 1 



Figure 3.1 The normalized beam intensityl/Io as a function of the radial distance 
p at different axial distances: (a) z = 0 (b) z = zo (c) z = 2zo [16] 


On the beam axis (p = 0) the intensity can be given by[15] 


/(0,z) = /c 




lf(z) 


I+(z/z,)^ 


(3.17) 
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It has maximum value Iq at z = 0 and drops gradually with increasing z, reaching half 
of its peak value at z = +Zo (figure 3.2). When|z| » z^, /(0,z) ?= /„ z,,'/^ so that the 

intensity decreases with the distance in accordance with an inverse-square law. The 
overall peak intensity I (0, 0) = Iq occurs at the beam center (z = 0, p = 0). 



Figure 3.2 The normalized beam intensity I/Io at points 
on the beam axis (p = 0) as a function of z [16]. 


Power 

The total optical power carried by the beam is the integral of the optical intensity over 
a transverse plane (at a distance z). 


P = ^l{p,z)27ipdp, P = ^loi^Wo) 


(3.18) 


Since beams are often described by their power P, it is useful to express /, in terms of 
Pas[15] 


IP 

l{p,z)= ,„ 2 :-Texp 
ttW {z) 


IfL 

W\z) 


(3.19) 


The ratio of the power carried within a circle of radius in the transverse plane at 
position z to the total power is 

2pl 


Pi) 

■ \l{p, z)2npd p = 1 - exp 


JV^(z) 


(3.20) 



Beam Radius 


Within any transverse plane, the beam intensity assumes its peak value on the beam 
axis, and drops by the factor ^2 = 0.135 at the radial distance p = W(z). Since 86% 

of the power is carried within a circle of radius fV(z) is regarded as the beam 
radius (also called the beam width). The rms width of the intensity distribution 

isa = — fV(z). The dependence of the beam radius on z is governed by 


W(z) = fV„ 



f \ 

2“ 


1 


1 + 

— 



^^0 ; 



1/2 


(3.21) 


It assumes its minimum value )To in the plane z = 0, called the beam waist. Thus Wo is 
the waist radius. The waist diameter 2lVo is called the spot size. The beam radius 
increases gradually with z, reaching ^/2^Fo, at z = zo, and continues increasing 
monotonically with z as shown in Figure 3.3. For z » zo, 

W 

W{z)«^z = e^z (3.22) 

^0 



Figure 3.3 Beam radius as a function of propagation distance.|161 

where 0o = Wo/zq. Using equation.(3.22) and equation.(3.1 1) we can also write 



Beam Divergence 

when z »zo, the beam radius increases approximately linearly with z, defining a cone 
with half-angle .About 86% of the beam power is confined within this cone. The 
angular divergence of the beam is therefore defined by the angle 


n 2W, 

The beam divergence is directly proportional 
1 and the beam-waist diameter 2fVo. 


(3.24) 

to the ratio between the wavelength 


Depth of Focus 

The axial distance within which the beam radius lies within a factor V2of its 
minimum value is known as the depth of focus. The depth of focus is twice the 
Rayleigh range. Since the beam has its minimum width at z = 0, as shown in 
figure 3.3, it achieves its best focus at the plane z = 0. In either direction, the beam 
gradually grows “out of focus.” 



Figure 3.4 The depth of focus of a Gaussian beam.116] 


Hence depth of focus is given by 
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Phase 

The dependence of phase of the Gaussian beam on the distance is given by 

= + (3.26) 

2R{z) 

On the beam axis (/? = 0) the phase is 


(p(Si,z) = kz~4{z) (3.27) 

Where kz phase of plan wave and ^(z) is phase retardation. 



Figure 3.5 Phase retardation of the Gaussian beam relative to a uniform 
planewave at points on the beam axis.[16] 


3.1.1 FWHM calculation by FRESNEL equation for the Gaussian Propagation 

To validate the FDBPM code we have compared the result of FDBPM with known 
result, so we have analyzed the propagation of Gaussian beam in free space. Using 
the Gaussian wave as an initial excitation with beam waist 3um we propagated upto 
400um. At 400um we have calculated the FWHM of the profile which is 5.02e-05um. 
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3.3 Error calculation in FWHM measured by FD BPM 


S.NO 

DELTAZ 

(meter) 

FWHM(FDBPM) 

(meter) 

FWHM(FRESNEL) 

(meter) 

ERROR 

(percentage) 

1 

l.OOOe-005 

4.158e-004 

4.158e-005 

1.717e+001 

2 

5.000e-006 

5.020e-004 

5.020e-005 

O.OOOe+000 

3 

2.000e-006 

5.020e-004 

5.020e-005 

O.OOOe+000 

4 

1 .OOOe-006 

5.020e-004 

5.020e-005 

O.OOOe+000 

5 

5.000e-007 

5.020e-004 

5.020e-005 

O.OOOe+000 


TableS.l Error in FWHM measured by FDBPM 


In the figure 3.8 the variation of the error in FWHM with the distance can be seen. 



Figure 3.8 Error in FWHM by Fresnel equation and FWHM by FDBPM 
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Figure 3.7 shows that FWHM measured by FDBPM is approaching the FWHM of the 
fresnel equation. Now form table 3.1 it was observed that with Az = the error was 

considerable less. So for further calculation I have taken Az = A . We have used the 
matlab to calculate the error which is accurately calculating floating point value upto 
32 total bits and 8 exponential bits. 


3.4 Mode Calculation Theory 

Since the Propagation constant in the z direction is given by exp(-//7z) so 
equation(2.14) chapter2 will become 


d^Ey 

dx^ 


= - CO^/XS]Ey 


A 



03 


(3.28) 


Figure 3.9 Slab Waveguide Structure 
Solution of the equation (3.28) can be expressed as 


= CicosA:ixX for -d<x<d 

(3.29) 

Ey(x) = Cl sin(A:i.xx) for -d <x< d 

(3.30) 

= co^/j£\ - 0^ =k\^ - 0^ 

(3.31) 

Ey(x) - D\ exp(-or 2 x) for x>d 

(3.32) 

a 2 ~ = -or usi = 0 -ki 

(3.33) 

Where C i and D i are constants. 



At the boundary both field are same so the field equation (3.32) will become 
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Ey{x) = Cicos(/:ixx)exp[-a2(j:-c?)] (3.34) 

At the boundary of two medium value of magnetic field Hu = Hu. The value of 
Hz can be calculated from equation(2.12) 

iklx ^ . y. iOC2 ^ y I *N 

Cism(A:i;cfi?) = C\cos(kixd) (3.35) 

COjU OJ/X 

Xdin{hxd) = — (3.36) 

ku 

Using sin(^u(i) as a solution again applying boundary condition Hu = H 2 z we get 

tan(kud) = -— (3.37) 

a2 

Using equation(3.31) , equation(3.33) we can write 

= A:o^[ni^ -n2‘]-^i;c^ (3.38) 

Normalized frequency variable V is given by 

= [m^-n2^]ko^d^ ( 3 . 39 ) 

Using equation (3.38) and equation (3.39) becomes 

=a2^d^+ku^d^ (3.40) 

We can get the solution of equation(3.36) and equation(3.37) by using the 
equation(3.40) . The number of solution of equation(3.36) and equation(3.37) have 
same number of modes will propagating in the slab waveguide. 

3.5 Verification of result for Waveguide 

Assuming the following data for the slab waveguide mode verification. 

Width of core =6um 
Core Refractive Index= 1.500 
Cladding Refractive Index=l .495 
Wavelength = lum. 
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Using Mathematical analysis with equation (3.36), equation (3.37) and 
equation (3.40) we get only one solution. Which shows that only one mode is 
propagating inside the slab waveguide as shown in figure 3.10. 



X A?(i$ ( 33 ym 


Figure 3.10 Amplitude Profile of Slab Waveguide structure by calculation 



XAxis( .33um ) 


Figure 3.11 Amplitude Profile of Slab Waveguide structure by FDBPM 
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Using the numerical method also we get the field profile as shown in the figure 3.1 1 
Which shows that we are getting the single mode propagation inside the slab 
waveguide. 

From the above analysis it can be concluded that FDBPM is suitable for waveguide 
structures. Which will be used for further analysis of the perturbation in waveguide. 
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Chapter 4 

RESULTS AND DISCUSSION 


During the first stage of the work, a FD BPM code considering the transversal 
electric field component is made and validated by free space Gaussian beam 
propagation. In the second stage, FDBPM is used for analysis of slab waveguide 
structure. Single mode and multimode analyses of slab waveguide structure are 
discussed in this section. Then, by FDBPM code, perturbed slab waveguide is also 
analyzed. The effect of different refractive indexes of the surrounding medium is also 
tested. 

4.1 Free Space Propagation 

Using FDBPM code we can see the intensity profile after launching the Gaussian 
excitation as an input for free space propagation. 

4.1.1 FDBPM without Boundary Condition 

As we analysis the physical structure numerically we can not analyze the structure for 
infinite width. We have to analyze the structure within the certain boundary. As we go 
out from the computational window their sudden change in the window so that due 
this change some reflection starts coming back into the computational window. 
Figure (4.1) shows the intensity profile using FDBPM without applying the 
boundary condition 
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ii 



Figure (4.1) Intensity profile without boundary condition 


4.1.2 FDBPM with Boundary Condition 

To remove the reflections coming from the boundary we apply the boundary 
condition. The following figure (4.2) shows the intensity profile using FDBPM with 
Transparent boundary condition 



Figure (4.2) Intensity profile with boundary condition 
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X Axis 


Figure (4.3) Intensity profile at 400um 

From figure(4.1) and figure( 4.2) we can see the significance of the boundary 
condition applied to the computational window. Figure (4.3) illustrates the 
interference pattern formed by the reflections in the presence of the computational 
boundaries and effect of TBC at 400um. 


4.2 Singlemode and Multimode Waveguide Structure Analyses 

Surroundine ^ 

Clad( 1.495) ^ 

Cored .500f 

37pm Clad(1.495) 



lOum Surrounding 



Figure (4.4) Slab Waveguide Structure 
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From section (3.5) chapter 3 we can see that existing code can be applied for slab 
waveguide structure. Using slab waveguide structure (figure (4.4)) the intensity 
profile of the single mode propagation was calculated, as illustrated in figure (4.5) . 



Figure 4.5 Single Mode Intensity profile 



Figure 4.6 Half of the Intensity Profile at Output 
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The effect of change in core width on the intensity profile of the structure was also 
analyzed. We varied the core width from 6/^mto M/um. We can see from the 
intensity profiles shown in figure (4.5) and figure (4.7) that single mode profile is 
converted into the multimode profile which is the effect of increase in the core width. 
It can inferred from this result that increase in the core width results the increase in 
number of mode. It was also observed that the amplitude of the intensity profile 
repeats in the periodic interval (figure (4.7)). 



Figure (4.7) Top View of the intensity profile using multimode structure 


Figure 4.8 illustrates the intensity profile for single mode and multimode structures 
while keeping surrounding medium as air at distance of 1 cm. 
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Figure (4.8) Intensity Profile of single mode and multimode structure 


4.3 Power Calculation 

In this thesis research we have tried to implement the refractive index sensors. This 
sensors works sense different refractive index accordingly change in output power. As 
we move along Z direction we get the intensity profile which is stabilized after 200um 
of propagation from the start of excitation. 

We can get the power P form the intensity profile I which is defined between interval 

X| and X 2 is defines as P = £ /c£c . Here we get the power form the intensity profile 

by numerical integration method. "Now a days there are two famous techniques used 
for numerical integration. 
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4.3.1 Trapezoidal Method 

The Trapezoidal Method is a numerical method for estimating the value of an integral 

£ f{x)dx or estimating the area under a curve generated by a function f(x). A linear 

function, which is the simplest polynomial approximation, is the basis of the 
Trapezoidal Method. As the name suggests, the Trapezoidal Method is based on 
taking the area of trapezoids after a given function f{x) has been approximated by a 
linear function 


Figure 4.9 Trapezoid area calculation for two points 

The area of the trapezoid is just the area of the rectangle plus the area of the triangle. 
That means our approximation is 

£ f{x)dx = (h - a)f(a) + 1 (Z» - a)[f{b) - f{a)] (4.1) 

^{h- a)[f{a) + i m - i /(«)] (4.2) 

= \{b-a)[m + fm (4.3) 

we cannot expect this to be a good approximation. However, we can break the region 

[a, 6] into many smaller pieces and apply the approximation on each piece. On the 

/ 

smaller pieces, the graph looks more and more like a straight line so the 
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approximation should improve. Let's choose some positive integer n and break the 
interval [a,b\ into n equal pieces. The width of each piece is /? = - — - . 



Figure 4.10 Trapezoid area calculation for more points 


We will label the points defined by the sub-intervals by Xj and call yi= ^ f{xi)clx . If 
we the approximate the area under the graph by the area of the trapezoids, we have 


f{x)dx ^Tn = ^{yo + y\) + ^{y\+yl) + ....^{yn-^ + yn) 


(4.4) 


= -{y^ + y\ + y\-\-y7 + yi+yi + + >] 


(4.5) 


- ”[yo "t" 2_yi + 2_y2 + 2y3 + + '2^yn - 1 + yn\ 


(4.6) 


4.3.2 Simpson’s Rule 

Simpson's rule is a approximating the integral of a function / using quadratic 
polynomials (i.e., parabolic arcs instead of the straight line segments used in the 
trapezoidal rule). Simpson's rule can be derived by integrating a third-order Lagrange 
interpolating polynomial fit to the function at three equally spaced points. In 
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particular, let the function / be tabulated at points xo, X|, and X2 equally spaced by 
distance h, and denote fn = f{xn) . Then Simpson's rule states that 


J pxi /•r(i+2/i 

f{x)dx= f{x)dx ( 4 . 7 ) 

xo J.to 

«^^(/» + 4/i + /2) (4.8) 

Since it uses quadratic polynomials to approximate functions, Simpson's rule actually 
gives exact results when approximating integrals of polynomials up to cubic degree. 

In exact form, 

f{x)dx = + 4/i+ /2) + -^ 


+ (4.9) 

x\ 

= 1;,(/o + 4/i + /2) + 7?« . (4.10) 

where the remainder term can be written as 

= — //'/<"’ (x*) (4.11) 

90 

with X* being some value of x in the interval [xi,X2] . 

An extended version of the rule can be written fory(jc) tabulated at Xo, X|, ...,X2n as 
f fix)dx = i h[fo + 4(/l + /j + + f2n-\) 

J.xn 3 

+2(/2 + /4+ + f2»-2) + f2n]-Rn (4.12) 

where the remainder term is 

Rn = — h ' f *\ x *) (4.13) 

90 

with X* being some value ofx in the interval [xi,X2k]. 
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4.4 Refractive Index Sensing with Evanescent Field 

Optical fiber is partially uncladded and is surrounded by a medium whose refractive 
index is to be measured. If optical fiber is partially uncladded but removal of cladding 
is not 100%, some cladding is remains in this case one part of this evanescent wave is 
transform into a propagating wave. By putting a third medium close to the 

interface, the evanescent wave is frustrated and a wave propagates into the third 
medium. This process is called frustrated total reflection [17]. It is given in 
figure 4. 1 1 . In Figure 4.11 d is the thickness of the cladding. 



Figure 4.1 1 Total electric field of the frustrated total reflection process in TE- 
mode. On the right, the profile in the z-direction is plotted. |17] 


4.4.1 Mathematical Analysis of Evanescent Wave Sensor 

Evanescent wave sensor works [18,19]. If a portion of the cladding is stripped away 
and a medium of different refractive index is placed in contact with the fiber 

core, and we find that for total internal reflection, the angle at which the ray entering 
at boundary between core and measured refractive index, must be greater than new 
critical angle: 


sin^; = «„,/«, 


(4.14) 





49 


Different guiding characteristics will apply and the new NA will be given by: 

cos^2 (4.15) 

For the case where core refractive index greater than measured refractive index, if the 
refractive index of the medium surrounding the exposed fiber core increases so that it 
approaches the refractive index of the core then the value of new critical angle 
increases and therefore new NA will decrease. This results in some of the light being 
lost from the fiber into the surroundings. 

When liquid is absorbing; the transmitted light through absorbing medium is given 
by; 

/ = /„ exp(-mZ.) (4.16) 

Where L is the length of absorption region and a is the absorption coefficient of the 
absorbing medium. The fraction of evanescent power, which can be, interacts with the 
absorbing medium is the quantity closely related to the fraction r, of the total guided 
power that resides in the cladded region, i.e. 

(4.17) 

Substantial values (>50%) for r iu single mode fibers and for multimode fiber the 
average fractional power in the cladding is very low (~1%). 


4.4.2 Sensitivity of the Sensor 

The sensitivity s of a function f(x) to parameter changes can be represented as follows 


[ 20 , 21 ]: 


^ Af{x)lfix) 
Aa,/a^ 


(4.18) 


Where a the parameter of interest are core diameter and cladding refractive index. 
The sensitivity s of V to changes in n 2 is given by 
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I' 



In^ dV 
V a«2 


2/72 ^ ^2 

^ A V«f -«2 


(4.19) 


From depth of penetration, the sensitivity s of as a result of change in n 2 is given 


by 




= ■ 

"2 


{Sin^e-Cy,)) 


(4.20) 


The negative sign in equation (4.19) indicates that an increase in the refractive index, 
n 2 , results in a decrease in the number of modes to be propagated into a fiber. The 
positive sign in (4.20) indicates that the penetration depth increases with cladding 
refractive index, a reduction in the electric field within the core. 


4.4.3 Factors affecting Sensitivity 

The various factors affecting the intensity of light lost are wavelength of the light, 
temperature, length of the fiber (cladded and un-cladded), shape of the fiber (micro 
bending) in the sensing uncladed region, effect of cladding thickness etc. In this 
experiment wavelength of light if fixed .The effect of temperature must be minimized 
to increase accuracy. 


4.5 Analysis of Perturbation in Slab Waveguide 

In this section first, the different percentage of cross sectional cut in some portion of 
the cladding is analyzed. Then its effect on Intensity profile and power at the output is 
observed. Effect of different refractive index is also analyzed at the output. 
Figure (4.12) and figure (4.13) shows the perturbation in the slab structure. The 
structure (figure(4.5)) is used as the reference for the power measurement and 
Intensity Profile calculation (figure(4.5)) at the output end. 



Refractive Index 
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Figure 4.12 Slab waveguide Structure with Cross Sectional Cut( 3d view) 
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Figure 4.13 Slab Waveguide Structure with Cross Sectional Cut (top view) 
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4.5.1 Effect of Cladding Removal 

Figure (4.14) illustrates the intensity profile during the propagation. The small 
increment in the cross sectional percentage cut is causing small decrement in the total 
output power of core and cladding (at the end). It is shown in the table (4.1) and 
figure (4.15). This decrement in the output power is due to leakage of more power 
coupled to the surrounding medium. 



Figure 4.14 Intensity Profile for 70 Percentage of cut in cladding 


S.NO 

DISTANCE 

(meter) 

CLADDING REMOVAL 

(In Percentage) 

POWER 

(watt) 

1 

3.000e-002 

00 

3.7599e-006 

2 

3.000e-002 

70 

3.7426e-006 

3 

3.000e-002 

80 

3.7329e-006 

4 


90 

3.726e-006 

5 

3.000e-002 

100 

3.4967e-006 

6 

3.000e-002 

100 percent cladding and 5 
percent core 

2.741 5e-006 


Table 4.1 Change in Output Power for Different percentage of cladding 
Removal(Surrounding medium Refractive Index=1.00) 











Figure 4.15 Intensity Profile for different Percentage of cut in cladding 


4.5.2 Effect of Refractive Index on Output Power 

Now we will discuss the effect of refractive Index on the output Power. As shown in 
the Figure (4.16) as we increase the refractive index of surrounding medium to the 
value up to cladding refractive index, the output power decreases because of decrease 
in the total internal reflection angle an so more power will go to the surrounding 


medium. 






Output Power (in log scale) 
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Figure 4.16 Effect of Refractive index (less than core) on the output power for 
different percentage of cladding removed structure 
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Figure 4.17 Effect of Refractive index on the output power for different 
percentage of cladding removed structure 
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Numerical Results shown in figure 4.16 also matches with practical results [22]. 
Now we will see the effect of refractive index more than the value of core refractive 
index (1.5) which is shown in Figure 4.17. From figure 4.19 we can see the combined 
effect of cross sectional cut and change in the refractive index. In each case of 
increase in the cut of cladding as well as the increase in refractive index, the output 
power becomes insensitive to the increase in the refractive index . 

4.6 Optical Sensor 

The new structure figure 4.18 is proposed here. In this structure two identical flat 
fibers are taken. Each flat fiber is allowing only single mode .They are separated 5um 
from each other. The 100 percent cross sectional cladding is removed upto 1cm from 
1cm of start of launching the power in fiber. Figure 4.23 shows that as we increase the 
refractive index of surrounding material the power is coupling to other waveguide . It 
also tells that when the refractive index of surrounding material is approaching the 
value of core refractive index (1.5) then power is more coupling in the surrounding 
medium and other waveguide. Figure 4.19 and figure 4.20 shows the intensity profile 
at refractive index of surrounding is less than 1.5 and Figure 4.21 and figure 4.22 
shows the intensity profile at refractive index of surrounding is 1.5. The sensitivity of 
this fiber with refractive index is shown in figure 4.23 . so Figure 4.23 can be used for 
measurement of optical properties of surrounding material like refractive index. 












Intensity Profile , , Intensity ftofile (in log) 


v/ / 


Figure 4.19 Intensity Profile for Surrounding Refractive Index at 1.44 


Figure 4.20 Intensity Profile for Surrounding Refractive Index at 1.44 



Figure 4.21 Intensity Profile for Surrounding Refractive Index at 1.50 
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Figure 4.22 Intensity Profile for Surrounding Refractive Index at 1.50 



Figure 4.23 Effect of Refractive Index on Output Power 
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4.7 Effect of Strip 

As fiber sensor is placed in liquid. Then liquid effects the power following in the 
fiber. Now it is assumed that liquid is not reacting with fiber so that no chemical 
change is coming in fiber or in the surrounding liquid. But if we put this sensor for 
more time in the liquid then small strip of layer is formed besides the fiber. This small 
strip affects the'power flowing in the fiber. Here 1 have analyzed the effect of strip on 
the output power. 1 have taken strip of width lum on both side of core at 100 percent 
clad removed section. Water is taken as surrounding material. Without any strip we 
are getting the output power 3.54t4e-006 watt. Now by using strip we can see from 
figure 4.25 that the output power is greater when refractive index of strip is near the 
refractive index of core. Table 4.2 and Figure 4.25 gives the output power at different 
refractive index of strip. Figure 4.24 is shows the structure used for analysis of strip 



X Axis 

Figure 4.24 Structure with lum strip on both side of core 


7. Strio 


Output Power 
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yio'® Effect Of strip 



Figure 4.25 Effect of Strip on Output Power with Surrounding Refractive Index 


Rl 

(Surr. Medium) 

Output Power(w) 
(with Strip) 

Output Power(w) 

(No strip) 

1.000 

3.585e006 


1.330 

3.492e-006 


1.400 

3.513e-006 


1.450 

3.695e-006 

3.5414e-006 

1.500 

3.678e-006 


1.550 

3.674e-006 


1.600 

3.433e-006 



Table 4.2 Effect of strip Refractive Index on Output Power 
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4.8 Effect of Roughness in cladding Removal 

Uptill now we have seen the effect of pure refractive index of surrounding material 
on the propagation of optical signal. We have also seen the effect of cladding 
removal. In real sensors when cladding is etched by enchants then surface is not 
perfectly flat etch. This roughness may effect the propagation of wave. In this section 
effect of roughness is analyzed. I have taken the sinusoidal perturbation at the core 
surface. The period of grating is equal to wavelength of propagation of optical signal. 
The structure used for the analyses of rough surface is shown in figure 4.26. 



Figure 4.26 Waveguide structure with roughness due to cladding removal 
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Figure4.27 Effect of Roughness and Refractive Index on Output Power 


RI 

(Surr. Medium) 

Output Power(w) 
(with Grating) 

Output Power(w) 
(with flat surface) 

1.000 

3.760e-006 

3.497e-006 

1.330 

3.760e-006 

3.541 e-006 

1.400 

3.760e-006 

3.582e-006 

1.450 

3.760e-006 

3.599e-006 

1.500 

3.723e-006 

2.287e-007 

1.550 

3.724e-O06 

9.685e-010 

1.600 

3.725e-006 

1.847e-009 


Table 4.3 Effect of Roughness and Refractive Index on Output Power 


Figure 4.27 shows the output power variation with roughness with different 
Refractive Index. Table 4.3 compares the output power with grating and without 
gratings for different refractive index of surrounding material. Table 4.3 shows that 
roughness at cladding is working as a grating which is trying to reduce the power 
leakage at the surrounding medium. This effect is clearly seen when the refractive 
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index of surrounding material ( table 4.3) is increased above core refractive index. 
Here we are using fiber structure for single mode Propagation. Which signifies that 
we are using only one wavelength for propagation. Roughness is working as grating 
which is reflecting leaking power to the core. Using such type of fiber with grating we 
can reduce the width of fiber with removing extra cladding and surrounding medium. 
Such Type of fiber can help to reduce the area of integrated chip 
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CHAPTERS 

CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 


The aim of this thesis is to develop the numerical code using Finite Difference Beam 
Propagation Method. In this work a detailed analysis of Numerical Method FDBPM is 
presented. FDBPM is applied to a perturbed slab waveguide structure and effect of 
cladding removal and effect of refractive index is analyzed. The results are concluded as 
follows: 

• An numerical method using FDBPM is developed for analysis of perturbed slab 
stmctures. Numerical code is also verified for free space Gaussian propagation and 
single mode slab waveguide structure 

In the Perturbed slab waveguide we get the following results 

• As the cross sectional cut percentage of the cladding increases the total output power 
of the core and the cladding is decreases. 

• On increasing the refractive index of the surroimding Medium to the value up to the 
cladding refractive index the output power decreases. 

• The combined effect of cross sectional cut of the cladding and change in the 
refractive index was observed with the increase of the cut of the cladding the output 
power become insensitive to the increase of the surrounding Medium refractive index 



The effect of strip besides the core is also observed as the refractive index of strip 
was near the core then output power is more than the without strip power. 


5.1 SUGGESTIONS FOR FUTURE WORKS 

Accuracy of FDBPM can be increased by Appling the Wide angle Approximation 
and second order Transparent Boundary Condition. 

The result obtained by using this method can be validated with the experimental 
data. 

In our work we consider only propagation of Gaussian beam, through a single slab 
waveguide structure. This work can be extended to array of waveguide structure 
analysis. 

In our estimation we assumed for the propagation of Gaussian beam has slowly 
varying envelop. It is helping to take large step size. It will be help full for saving 
the memory space as well as time of computation. It can also be extended for fast 
varying envelop to get more accuracy. 
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